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Abstract. We discuss steady-state transonic outflows ob- 
tained by direct numerical solution of the hydrodynamic 
and magnetohydrodynamic equations. We make use of the 
Versatile Advection Code, a software package for solving 
systems of (hyperbolic) partial differential equations. We 
proceed stepwise from a spherically symmetric, isother- 
mal, unmagnetized, non-rotating Parker wind to arrive 
at axisymmetric, polytropic, magnetized, rotating mod- 
els. These represent 2D generalisations of the analytical 
ID Weber-Davis wind solution, which we obtain in the 
process. Axisymmetric wind solutions containing both a 
'wind' and a 'dead' zone are presented. 

Since we are solving for steady-state solutions, we effi- 
ciently exploit fully implicit time stepping. The method al- 
lows us to model thermally and/or magneto-centrifugally 
driven stellar outflows. We particularly emphasize the 
boundary conditions imposed at the stellar surface. For 
these axisymmetric, steady-state solutions, we can use the 
knowledge of the flux functions to verify the physical cor- 
rectness of the numerical solutions. 
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1. Introduction 



solar wind 



Observational and theoretical research on stellar winds 
and astrophysical jets has evolved rapidly. For our own sun 
and its associated solar wind, the current understanding 
necessitates the combined study of solar wind acceleration 
and coronal heating in time-dependent modeling (Holzer 
& Leer 1997). At the same time, Holzer and Leer right- 
fully stressed that it remains useful to emphasize on early 
studies of wind acceleration. This, we find, is especially 
true for numerical modeling of stellar winds. With the ul- 
timate goal of time-dependent heating/wind modeling in 
mind, we here address the simpler question on how to ac- 
curately model ID and 2D steady-state winds by the nu- 
merical solution of the polytropic magnetohydrodynamic 
(MHD) equations. 

Send offprint requests to: R. Keppens 



Since much of the ID solutions we obtain is known 
from the outset, we can verify our results precisely. Indeed, 
ID Parker Ql95S| ) and Weber-Davis (Weber & Davis |1967D 



wind solutions can be checked to agree with the analytic 
description. In their axisymmetric 2D extensions, vari- 
ous physical quantities must be conserved along poloidal 



streamlines. Sakurai (1985, 1990) presented such 2D gen- 



eralizations of the magnetized Weber-Davis wind, using 
a method designed from these conservation laws. With 
modern numerical schemes, we can recover and extend 
his solutions to allow for the self-consistent modeling of 
'dead' and 'wind' zones, as in the solar wind. Our steady- 
state solutions can be checked to conserve quantities along 
streamlines a posteriori. 

Of crucial importance is the choice of boundary con- 
ditions used in the simulations. Since the governing equa- 
tions for steady-state, transonic MHD flows are of mixed- 
type, their character can change from elliptic to hyper- 
bolic at a priori undetermined internal critical surfaces. 
Causality arguments have been used to discuss which and 
how many boundary conditions must be prescribed (Bo- 
govalov 1997). Our choice of boundary conditions used at 
the stellar surface is therefore discussed in detail. 

All solutions presented are obtained with a 
single software package, t he Ve rsatil e Advec- 
tion Code (VAC, see T oth |l996t |l997| and also 
tittp: //www. phys.uu.nl/~toth). Although we only 
present steady-state transonic outflows in spherical and 
axisymmetry, VAC is developed for handling hydrody- 
namic (HD) and magnetohydrodynamic (MHD) one-, 
two-, or three-dimensional, steady-state or time-depen- 
dent problems in astrophysics. It is therefore capable of 
achieving our ultimate goal of time-dependent 3D wind 
modeling. The insight gained in this study of steady-state 
polytropic flows will then be very useful. 

In Sect. ||, we list the equations and discuss the Versa- 
tile Advection Code for solving them in Sect. ||. Our cal- 
culations are presented in Sects. [|, [| and |6[ Conclusions 
are given in Sect. [7| The approach taken is a gradual one, 
where for instance our ID solutions are used to construct 
initial conditions for their 2D extension. We will therefore 
model, in increasing order of complexity: (i) isothermal, 
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spherically symmetric Parker winds; (ii) polytropic, spher- 
ically symmetric Parker winds; (iii) polytropic, rotating 
Parker winds for the equatorial plane; (iv) Weber-Davis 
magnetized, polytropic, rotating winds for the equatorial 
plane; and finally axisymmetric, polytropic, rotating 2D 
winds, both (v) unmagnetized and magnetized, without 
(vi) and with (vii) a 'dead' zone. 

2. Equations 

We solve the HD and MHD equations expressed in the con- 
servative variables density p, momentum vector pv, and 
magnetic field B. These are given by 



! + V>v) = 0, 



dt 

dB 

~dt 



V • [pvv + p tot l - BB] = pg, 



V • (vB - Bv) = 0. 



(1) 



(2) 



(•3) 



We introduced ptot = p+\B 2 as the total pressure, I as the 
identity tensor, g as the external gravitational field, and 
exploited magnetic units such that the magnetic perme- 
ability is unity. We drop the energy equation and assume 
a polytropic relation connecting the thermal pressure p 
and the density p. For a polytropic index 7, we thus as- 
sume p ~ p 1 . Hence, we do not address the heat deposi- 
tion in the corona. Although we solve the time-dependent 
equations as given above, we will only present steady-state 
d/dt = solutions of Eqs. For stellar wind cal- 

culations, we consider a spherically symmetric external 
gravitational field g = -GM,/r 2 e n where G is the gravi- 
tational constant, M» is the stellar mass, r is the distance 
to the stellar center, and e r indicates the radial unit vec- 
tor. 



3. Versatile Advection Code 

In this section, we discuss the software and numerical 
method used. The physics results are described from 
Sect. ^ onw ards. The Versatile Advection Code (VAC, see 
Toth 1996 , 1997 ) is a general purpose software package 
for solving a conservative system of hyperbolic partial dif- 
ferential equations with additional non-hyperbolic source 
terms, such as the MHD equations. VAC runs on PC's, 
on a variety of workstations, on vector platforms, and we 
can also run in parallel on a cluster of workstations, and 
on distributed memory architectures like the Cray T3D 



and T3E, a nd the IBM SP (Keppens & Toth |1998| , Toth 
& Keppens 1998). The code is written in the dimension 
independent LASY syntax (Toth [1997 ), so it can be used 
as a convenient tool to handle HD and MHD one-, two-, 
or three-dimensional problems in astrophysics and labo- 
ratory plasma physics. The dimensionality of the problem 



and the actual set of equations to solve are easily selected 
in a preprocessing step. 

VAC uses a structured finite volume grid and offers 
a choice of conservative, second order accurate, shock- 
capturing, spatial and temporal discretization schemes. 
The spatial discretizations include two Flux Corrected 
Transport variants and four Total Variation Diminishing 



(TVD) schemes (Toth & Odstrcil |1996J). Temporal dis- 
cretization can be explicit, semi-implicit, or fully implicit. 



It was recently demonstrated (Keppens et al. 1998, Toth 



et al. 1998) how the implicit approach can be used very ef- 



ficiently, for steady-state and time-accurate problems pos- 
sibly containing discontinuities. Here, we expect smooth 
solutions to the steady-state HD and MHD equations, so 
one can greatly benefit computationally from fully implicit 
time stepping. 

In this paper, we solve the polytropic HD (B = 0) 
Eqs. (|) and (§) in Sects, fy] and & 

Magnetohydrody- 

namic equations are solved in Sects. 12 and ||. We solve 
one-dimensional problems in Sect. || and two-dimensional 
problems in Sects. [5] and |[ In practice, this means that the 
stellar winds we model are solutions of the equations under 
the additional assumption of a prescribed symmetry in the 
ignored directions. One-dimensional problems assume a 
spherical symmetry, while 2D solutions assume d/dip = 0. 
Here, ip denotes the angle in a cylindrical (R, <p, z) co- 
ordinate system centered on the star with its polar and 
rotation axis as z-axis. 

Since we are interested in steady-state solutions, we 
use fully implicit time stepping as detailed and demon- 
strated in Keppens et al. ( |1998| ) and Toth et al. ( |1998[ ). 
The linear systems arising in the linearized fully implicit 
backward Euler scheme are solved using a direct block 
tridiagonal solver for the ID problems and using a precon- 
ditioned Stabiliz ed Bi -Conjugate Gradient iterative solver 
(van der Vorst |l992| ) for the 2D cases. The Modified 



Block Incomplete LU preconditioner is described in van 
der Ploeg et al. (|l997| ). We consistently used the TVD 



Lax-Friedrich (T VDLF ) spatial discretization (Yee |1989 
Toth & Odstrcil 1996| ) using Woodward limiting (Collcla 



Woodward 1984 ). We typically took Courant numbers 
C = 0(100). For all ID solutions and for the 2D hydro- 
dynamic solutions, the steady-state is reached when the 
relative change in the conservative variables from one time 
level to the next drops below 10 -8 . We use a normalized 



measure defined by (Toth et al. 199S 



A 2 U = 



\ 



iV var 

E 



S K rid(^« 



n+1 



us) 2 



n\2 



(4) 



where N var is the number of conserved variables Ui, and 
the superscripts indicate the pseudo time level t n . Since 
we are solving for smooth solutions, the numerical schemes 
can easily achieve such accuracy in the steady-state solu- 
tions. 
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For the axisymmetric 2D MHD solutions, the way to 
ensure a zero divergence of the magnetic field is non- 
trivial. Unless one uses a scheme which keeps V • B = ex- 
actly in some discretization, like the constrained transport 



method (Evans & Hawley 1988), corrective action needs to 
be taken during the time integration. This involves either 
including corrective source terms in the equations which 
are proportionate to the numerically generated divergence 



(Powell 1994, their use f or TVDLF was first advocated 



by Toth & Odstrcil 1996 ), o r mak ing use of a projection 
scheme (Brackbill & Barnes 198C ) which involves the so- 
lution of a Poisson equation. These two approaches can 
even be combined and such combination is also beneficial 



for fully implicit schemes (Toth et al. g_998J). For the 2D 
axisymmetric MHD wind solutions presented in Sect. ^, 
it proved difficult to ensure a divergence free solution in 
a fully implicit manner. Using explicit time stepping and 
employing the projection scheme before every time step, 
we could get steady-state solutions where ~ 0(1O -7 ) 
which have acceptable | VB |< 10~ 3 , although this time- 
marching method is computationally much more costly 
than the implicit approach. The use of Powell source terms 
alone proved inadequate for these wind solutions, as fairly 
large errors are then advected into the whole solution do- 
main from the stellar boundary outwards. 

VAC makes use of two layers of ghost cells surrounding 
the physical domain to implement boundary conditions. 
A symmetry condition at a boundary is then imposed by 
mirroring the calculated values in the two physical cell 
layers adjacent to the boundary in these ghost cells. The 
boundary conditions imposed at the stellar surface are ex- 
tremely important. We will emphasize them for all cases 
considered. 

4. ID polytropic winds 

4-1- Parker winds 

Our starting point is the well-known analytic Parker 
( |1958| ) solution for a spherically symmetric, isothermal 
(7 = 1) outflow from a star of mass Af» and radius 
r„. Given the magnitude of the escape speed v esc = 



-y/2GAf»/r*, one can construct a unique 'wind' solution 
which starts subsonically at the stellar surface and accel- 
erates monotonically to supersonic speeds. This solution is 
transonic at the critical position r s — v 2 sc r » / '(4c^) , where 
c l* = p/p i s the constant isothermal sound speed. Since 
we know the position of the critical point r s , we can eas- 
ily determine the flow profile v r (r) and the corresponding 
density profile p(r). The radial velocity is obtained from 
the iterative solution of the transcendental equation 



\ 



r 



3 + 21n 



(5) 



The density profile results from the integrated mass con- 
servation equation. Since the radial velocity reaches a con- 



stant supersonic value asymptotically, the corresponding 
density vanishes at infinity as 1/r 2 . 

Choosing units such that r* = 1, p* = p(r*) = 1 with 
c s * = 1, we initialize a ID spherically symmetric, poly- 
tropic (with 7 > 1) hydrodynamic outflow with this ana- 
lytic isothermal Parker wind solution on a non-uniform 
mesh ranging through r £ [l,50]r». We use 1000 grid 
points and exploit a grid accumulation at the stellar sur- 
face, where the acceleration due to the pressure gradient 
is largest. In the ghost cells used to impose boundary con- 
ditions at the stellar surface, we fix the value of the base 
density to unity, and extrapolate the radial momentum 
continuously from its calculated value in the first grid cell. 
At r = 50r*, we extrapolate both density and radial mo- 
mentum continuously into the ghost cells. We then use a 
fully implicit time integration to arrive at the correspond- 
ing steady-state, spherically symmetric, polytropic Parker 
wind solution. The obtained solution p(r) 1 v r (r) can be 
verified to have a constant mass flux pr 2 v r as a function 
of radius and energy integral 



2 

s* ' 



E = (v 2 /2 + c 2/( 7 - 1) - GAU/r) /c; 

where c 2 (r) = p 7_1 . Also, determining the sonic point r s 
where v r (r s ) = c s (r s ) and the base radial velocity v r *, the 
solution can be checked to satisfy 



2(7 +1) 2(7 -1) 

Vesc \ 5-37 ( 2f£\ 5-37 
2c s * 



(0) 



Note how the isothermal 7 = 1 case is the only polytropic 
wind solution where the position of the critical point can 
be determined a priori. 

In practice, we increased the polytropic index gradu- 
ally from 7 = 1 through 1.05, 1.1, 1.12, 1.125 to 7 = 1.13, 
each time relaxing the obtained steady-state solution for 
one polytropic index to the unique transonic wind solu- 
tion for the next value. In the top panel of Fig. [j], we plot 
the radial variation of the Mach number M s — v r /c s for 
the isothermal 7 = 1 Parker wind with v esc = 3.3015c,;,,,, 
and for similar polytropic winds with 7 = 1.1 and 1.13. 
The vertical dashed lines indicate the agreement of the 
positions of the sonic points where M s — 1 with the cal- 
culated right hand side of Eq. (^J) . Note the outward shift 
of the sonic point with increasing polytropic index and the 
corresponding decrease of the asymptotic radial velocity. 

When we relax the restriction of spherical symmetry 
by allowing a rigid stellar rotation rate 57*, we can easily 
construct a solution for the equatorial plane only. Indeed, 
ignoring variations perpendicular to this plane, one simply 
adds a toroidal velocity profile where v v (r) = Q,*r*(r*/r) 
and then solves for p(r), v r (r), and v v (r). The boundary 
conditions on the toroidal momentum keep pv v fixed at 
the base and extrapolate it continuously at the end of 
the computational domain. A polytropic, rotating Parker 
solution for the equatorial regions is found by relaxation 
from a non-rotating wind with the same polytropic index 
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Fig. 1. Polytropic Parker winds with v esc = 3.3015c s *. 
Top panel: Mach number M s as a function of radial co- 
ordinate for spherically symmetric Parker winds for poly- 
tropic index 7=1 (isothermal), 7 = 1.1 and 7 = 1.13. 
Bottom panel: equatorial solutions for rotating 7 = 1.13 
polytropic Parker winds for various rotation parameters 
Critical points are indicated. 



7. In the bottom panel of Fig. [l], we show the Mach number 
M s (r) for Vesc = 3.3015c s * and 7 = 1.13 Parker winds 
where C = fi*r*/c s * equals ( — 0.0156, ( — 1 and ( = 1.9. 
The solution with ( = 0.0156 hardly differs from its non- 
rotating thermally driven counterpart shown in Fig. |l|, as 
expected. The additional centrifugal acceleration causes 
an increase in the base velocity v rSf and in the asymptotic 
radial velocity. Again, the solution can be verified to have 
a constant radial mass flux pr 2 v r , Bernoulli function 



E = (k 2 + v%]/2 + cJ/( 7 - 1) - GMjr) /c 2 s 



and constant specific angular momentum rv v . The posi- 
tions of the critical point (s) r s are now obtained from a 
generalization of Eq. (^) , namely from the solutions of 

2(7-1) 
7 + 1 



2 

v esc r * 



2c 2 



= 0.(7) 



This equation reduces to a second degree polynomial for 
a 7 = 1 isothermal, rotating, Parker wind so it is evident 
that rotation rates exist that introduce a second critical 
point. In Fig. [j], only the C, = 1.9 solution exhibits two 
critical points, shown as vertical dotted lines, within the 
domain. We determined the critical point(s) by solving 
Eq. (0) using the calculated base speed iv*. The close- 
up of the radial variation of M s for ( — 1.9 at the base 
reveals that this thermo-centrifugally driven wind passes 
the first critical point while being decelerated, then starts 
to accelerate and finally becomes supersonic at the second 
critical point. To correctly capture the dynamics close to 
the stellar surface it is clear that we need a high grid 
resolution, especially at the stellar surface. Indeed, the 
first critical point for the £ = 1.9 solution is situated at 
1.006rv 



.2. Weber- D 



avis win 



ds 



The magnetized Weber-Davis (WD) solution (Weber & 
Davis 1967) represents a valuable extension to the ro- 



tating, polytropic Parker wind solution for the equatorial 
plane. Again assuming that there is no variation perpen- 
dicular to this plane, one now needs to solve for an ad- 
ditional two magnetic field components B r {r) and B v {r). 
One is trivially obtained from the V • B = equation, 
namely B r — B r *r 2 /r 2 . The analytic treatment reveals 
that the magnetized polytropic wind solution has a to- 
tal of two critical points, namely the slow r s and the fast 
rf critical point. These are determined by the zeros of 
vf. - v 2 (c 2 + A 2 r + A 2 ,) + c 2 A 2 .. In between lies the Alfven 
point ta , defined as the radius at which the radial velocity 
v r equals the radial Alfven speed A r ~ B r / \fp. Since the 
equatorial fieldline is prescribed to be radial in the poloidal 
plane and the transfield force balance is not taken into ac- 
count, this Alfvenic transition is not a critical point in this 
model. 

In the fully implicit time stepping towards a steady- 
state WD wind for specific values of v esc — 3.3015c s *, 
7 = 1.13, ( = 0.0156, and for the base radial Alfven speed 
A r * = B r *j ^fpl = 3.69c s „, we initialize p{r), v r (r), and 
Vip{r) with the corresponding non- magnetic, polytropic ro- 
tating Parker solution. We fix B r to its known 1/r 2 de- 
pendence throughout the time evolution, and initialize B v 
to zero. The boundary conditions at r = 50r* extrapo- 
late all quantities we solve for continuously into the ghost 
cells. At the base, we keep the density fixed, the radial 
momentum and toroidal field component are extrapolated 
linearly from the first two calculated mesh points, while 
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the toroidal momentum pv v is coupled to the magnetic 
field ensuring 



10.00 F 



fLr* 



(8) 



This expresses the parallelism of the velocity and the mag- 
netic field in the frame rotating with the stellar angular 
velocity Q*. Using these initial and boundary conditions, 
we arrive at the unique WD wind solution for the given 
parameters. This magnetized polytropic wind solution for 
the equatorial plane is shown in Fig. |^. The solution agrees 
exactly with the analytic WD wind: we obtain five con- 
stants of motion, namely the mass flux pr 2 v r ~ 0.0139, 
the magnetic flux r 2 B r = 3.69 which is constant by con- 
struction, the validity of Eq. (|J) over the whole domain, 
the Bernoulli integral 

E = ([v 2 + v$\/2 + c 2/( 7 _ l) _ GM*/r - v v B v B r /pv r 
+B*/p) /<*, ~ 2.45, 

and the constant total specific angular momentum L = 
rv v — rB v B r j pv T ~ 13.36. The positions of the criti- 
cal points are r s = 7.4r* and r/ = 31.2r„«, while the 
Alfven point is at ta = 29.2r», as indicated in Fig. 0. 
This agrees with the values given in the Appendix to 
Keppens et al. ( 1995| ), where the same WD solution 
was calculated in a completely different fashion. Indeed, 
the WD solution for given values of 7, c s „, v esc /c s ^, C, 
and A„/c s „, can alternatively be calculated as a mini- 
mization problem in a six-dimensional space (see Belcher 
& MacGregor 1976) where one solves for the six un- 
knowns [v r *, v^, r s , v r (r s ), rf, v r (rf)]. This can be done 
using standard Newton-Raphson iteration provided one 
has an educated initial guess, but it can even be ob- 
tained by the use of a ge netic algorithm, as first demon- 
strated by Charbonneau ( 1995 ). The fact that the, for our 
method initially unknown, base velocities appear in the 
determining set of variables for these minimization meth- 
ods again indicates that a high base resolution is abso- 
lutely essential. Values for these six unknowns found from 
the solution in Fig. | are [0.01395,0.01541,7.4,0.6018, 
31.2, 1.1592] and these agree with the Newton-Raphson 
solution. 

The calculation of the WD polytropic wind by the 
stepwise relaxation from an isothermal Parker wind is thus 
an excellent test for the numerics, as every step (from 
isothermal to polytropic, from non-rotating to rotating, 
from Parker to WD) can be verified precisely to agree 
with the known solutions. It should be clear that we can 
construct WD wind solutions where the acceleration re- 
sults from the combined action of thermal, centrifugal, 
and magnetic forces. However, our interest is in the gen- 
eralization of these ID models out of the equatorial plane. 
We will again proceed in logical steps towards this goal. 
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Fig. 2. Weber-Davis wind solution for the equatorial 
plane. Top panel: radial variation of the radial Alfven 
speed A r , sound speed c s , and velocities v r and all 
normalized to the base sound speed c s *. Bottom panel: the 
corresponding poloidal Alfven Mach number Ma = v r /A r , 
and poloidal slow and fast Mach numbers, determining the 
positions of the critical points and the Alfven point. See 
text for parameters. 



5. Axisymmetric 2D polytropic HD winds 

To arrive at a crude model for the coronal expansion of 
a rigidly rotating star, we set forth to construct an ax- 
isymmetric, steady-state, polytropic wind solution valid 
throughout a poloidal cross-section. With the polar axis 
as rotation and symmetry axis, we need to generalize the 
rotating, polytropic Parker wind which succesfully mod- 
eled the equatorial regions. Whereas the Parker solution 
had at least one critical point, its 2D extension is expected 
to give rise to critical curves in the poloidal plane. The 
degree of rotation determines the deviation from perfect 
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circles arising in the non-rotating, spherically symmetric 
case. 

To initialize a 2D fully implicit time-stepping proce- 
dure to arrive at a steady-state wind, we use the ID 
Parker solution with identical escape speed v esc , poly- 
tropic index 7, and rotational parameter £. We use a 
spherical (r, 8) grid in the poloidal plane, where the grid 
spacing is equidistant in 9, but is accumulated at the 
base in the radial direction. We take a 300 x 20 grid and 
only model a quarter of a full poloidal cross-section. The 
density is initialized such that for all angles 8, the ra- 
dial variation equals the ID Parker wind appropriate for 
the equator. Writing the Parker solution as p p (r), v P (r), 
v P (r), we set p(r 7 8;t = 0) = p p {r), and similarly, we set 
v r (r,6;t = 0) = v P (r) and v v {r,6;t = 0) = v P (r)sm(8) 
so that it vanishes at the pole 6 = 0, while vg(t = 0) = 
everywhere. Since we now use a coarser radial resolu- 
tion, we interpolate the Parker solution linearly onto the 
new radial grid. Boundary conditions then impose sym- 
metry conditions at the pole (8 = 0) and the equator 
(8 = 7r/2). The radial coordinate still covers r £ [1, 50]r„, 
as in the ID calculations. Since the solutions are super- 
sonic at r = 50r„, the boundary conditions there merely 
extrapolate the density and all three momentum compo- 
nents linearly in the ghost cells. The stellar rotation en- 
ters as a boundary condition in the toroidal momentum 
component, which enforces v v = 0*7?*, where (R,z) are 
the cartesian coordinates in the poloidal (r, 8) plane. Note 
that the toroidal momentum may still change in the pro- 
cess, since we can no longer fix the density at the stellar 
surface to a ^-independent constant value. This is because 
in steady-state, the density profile should establish a gra- 
dient in the 8 direction to balance the component of the 
centrifugal force in that direction. In the purely radial di- 
rection, the inwards pointing gravity must be balanced by 
the combination of the pressure gradient and the radial 
component of the centrifugal force. We therefore extrapo- 
late the density linearly at the base. To enforce the total 
mass flux as in the equatorial Parker solution, we deter- 
mine the constant f mass = p p r 2 v P from the ID calcula- 
tion, and fix pv R = f mass R/r 3 and pv z = f m as S z/r 3 at 
the stellar surface for its 2D extension. 

An elementary analytic treatment for a 2D polytropic 
steady-state wind solution proceeds by noting that mass 
conservation is ensured when the poloidal momentum is 
derived from an arbitrary stream function x(i?, z) such 
that pv p = (1/R)e v x V%. It is then easily shown that 
the toroidal momentum equation is equivalent with the 
existence of a second arbitrary function L(x) — RVus, cor- 
responding to the conservation of specific angular momen- 
tum along a poloidal streamline. Similarly, energy conser- 
vation along a streamline introduces 

E(X) = [Vr + v 2 z + v%] /2 + pT-7( 7 - 1) - GAU/r. 

Across the poloidal streamlines all forces must balance 
out. 



We show streamlines and the contours of constant po- 
loidal Mach number M p = J [y\ + v 2 )/c 2 for two hydro- 
dynamic wind solutions for v esc = 3.3015c s *, 7 = 1.13, 
and with ( = 0.0156 (top panel) and £ = 0.3 (bot- 
tom panel) in Fig. 0. We restricted the plotting region 
to about 10rv For the imposed mass flux parameter, we 
used the values fmass — 0.01377 for the slow rotator and 
fmass = 0.01553 for the faster rotator, as found from the 
equatorial Parker solution for the same 7 and C- Note how 
for low rotation rates, the wind solution is almost spher- 
ically symmetric with nearly radial streamlines and cir- 
cular Mach curves. For higher rotation rates, the critical 
Mach curve where M p = 1 moves inwards at the equa- 
tor and outwards at the pole when compared to a non- 
rotating case. The streamlines show the equatorward de- 
flection when material is released from the stellar surface 
due to the centrifugal force. Such equatorward stream- 
line bending due to rotation is discusse d in d etail in the 
analytical study by Tsinganos & Sauty (1992). For these 
solutions, we can then verify that the specific angular mo- 
mentum L, as well as the total energy E, are conserved 
along the streamlines. 



6. Axisymmetric 2D polytropic MHD winds 

6.1. Magnetized winds 

To obtain an axisymmetric magnetized wind solution, we 
may simply add a purely radial magnetic field to a 2D 
HD wind solution and use this configuration as an ini- 
tial condition for an MHD calculation. Hence, we set 
B R (R,z;t = 0) = (iR/r 3 and B z (R,z;t = 0) = (3z/r 3 
while B v (t = 0) = 0. Such a monopolar field is rather 
unrealistic for a real star, but it is the most straightfor- 
ward way to include magnetic effects. The same type of 
field was used by Sakurai ( |1985| , |1990| ), which contained 
the first 2D generalization of the WD model. 

Boundary conditions at equator and pole are imposed 
by symmetry considerations. At r = 50r„, we extrapo- 
late all quantities linearly. Similarly, the base conditions 
extrapolate the density profile p and all magnetic field 
components, while the poloidal velocity components are 
set to ensure a prescribed mass flux. The stellar rotation 
rate and the coupling between the velocity and the mag- 
netic field enters in the boundary condition at the stellar 
surface where we demand 



B u 



vl/JBl + B 2 . 



(9) 



Specific attention is paid to ensuring the V • B = con- 
dition. As explained in Sect. |^, we now switch strategy 
and use explicit time stepping combined with a projec- 
tion scheme to obtain the steady state solution. 

We calculated the 2D extension of the WD wind cor- 
responding to v esc = 3.3015c s », 7 = 1.13, C = 0.0156 and 
(3 = 3.69. The mass flux is set to be f mass = 0.01377. 
We show in Fig. H the streamlines, and the positions of 



R. Keppens & J. P. Goedbloed: Numerical simulations of stellar winds 



7 



10H 




8 10 



10 £_-!— I 




Fig. 3. 2D Polytropic HD winds. We show streamlines and 
contours (dotted for values below unity) of the poloidal 
Mach number M p in the poloidal plane. For low (top 
panel) and high (bottom panel) rotation rates. 



the critical curves where the poloidal Alfven Mach num- 
ber and the poloidal slow and fast Mach numbers equal 
unity. The squared poloidal Alfven Mach number M V A is 
given by 



{M\f 



/(A 



Al 



(10) 



with Alfven speeds Aj = Bi/^Jp. The squared poloidal 
slow Mf and fast Mj Mach numbers are defined by 



(MfY 



2(4 



c « + Ap + Ay 



A 2 



4 c 2 A 2 



(11) 



A 2 

v 



Al 



Ac 2 s Al 



(12) 



At the pole, the fast Mach number coincides with the 
Alfven one since B v vanishes there and the parameters 
are such that A 2 > c 2 . Away from the pole, the toroidal 
field component does not vanish, so that Alfven and fast 
critical curves separate. Note how the equatorial solution 
strongly resembles the WD wind solution for the same pa- 
rameters shown in Fig. |^. The obtained wind solution is 
mostly thermally driven, like the solar wind. The rotation 
rate and magnetic field effects are minor an d an a lmost 
spherically symmetric wind results. Sakurai ( 1985 , 1990 ) 
demonstrated that for stronger fields, the magnetic force 
of the spiraling fieldlines deflect the outflow poleward. 
This magnetic pinching force can produce a polar collima- 
tion of the wind. These effects have also been addressed 
by analyt ical studies of self-similar outflows in Trussoni et 
al. ( |1997D . 

For these axisymmetric, steady-state MHD outflows, 
the solutions can be verified to obey the following conser- 
vation laws. Mass conservation is ensured when writing 
the poloidal momentum vector as pw p = (1/R)e v x V%, 
with the stream function x(-R, z). The zero divergence of 
the magnetic field yields, likewise, Bp = (l/i?)eL x V^, 
with ip the flux function. The poloidal part of the in- 
duction equation then leads to x("0)j provided that the 
toroidal component of the electric field E v vanishes. This 
can easily be checked from vrB z = v z Br, and the solution 
shown in Fig. ^ satisfies this equality to within 1%. This 
allows us to write x' = dx/dip — pv p /B p = pvr/Br = 
pv z /B z . A fair amount of algebra shows that the toroidal 
momentum and induction equation introduce two more 
flux functions, namely the specific angular momentum 
L(tp) = Rv v — RB v B p / pv p and a function related to the 
electric field fl(ip) = [v v — {v p /B p )B ip \ / R. The Bernoulli 
function derivable from the momentum equation can be 
written as 

E( X ) - [v i R + vt + v*]/2 + p'- i -/(>y-i)-GM./r 



vipBtpBp/ pv p 



Blip. 



Note how the constants of motion found in the WD solu- 
tion immediately generalize in this formalism (mass flux, 
magnetic flux, corotation as in Eq. (||), specific angular 
momentum L and Bernoulli function E). The hydrody- 
namic limit is found for zero magnetic field B = and 
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B z = a d 



(2z 2 - R 2 ) 
(R 2 + z 2 ) 5 / 2 ' 



(14) 



Fig. 4. Polytropic axisymmetric MHD wind. We show the 
streamlines and the positions of the critical surfaces where 
the poloidal Mach numbers equal unity. Parameters are as 
in the ID Weber-Davis wind shown in Fig. 0. 



vanishing electric field O = 0. Across the poloidal stream- 
lines, the momentum balance is governed by the gener- 
alized, mixed-type Grad-Shafranov equation. The numer- 
ical solutions we obtained indeed have parallel poloidal 
streamlines and poloidal fieldlines and conserve all these 
quantities along them. 

6.2. Winds containing a 'dead' zone 

The monopolar field configuration used above is unrealis- 
tic. However, it should be clear that our method easily gen- 
eralizes to bipolar stellar fields by appropriately changing 
the initial condition on the magnetic field. In fact, a star 
like our sun has open fieldlines at both poles and closed 
fieldlines around its equator. To obtain a steady-state stel- 
lar wind containing a 'wind' zone along the open fieldlines 
and a 'dead' zone about the stellar equator, we can sim- 
ply initialize the polar regions up to a desired polar angle 
Owtnd as above. The equatorial 'dead' zone is then initial- 
ized as follows: the density and the toroidal momentum 
component is taken from the 2D HD wind with the same 
rotational and polytropic parameters while the poloidal 
momentum components are set to zero. The initial mag- 
netic field configuration in the 'dead' zone is set to a dipole 
field which has 



B, 



3a d 



zR 



(R 2 +Z 2 ) 5 / 2 ' 



(13) 



The strength of the dipole is taken = (3/(2 cos(9 W i n d)) 
to keep the radial field component B r constant at 9 = 
9 wind- The initial B v component is again zero through- 
out. In summary, we now have the following set of pa- 
rameters used in the simulation: the escape speed v esc , 
the polytropic index 7, the rotational parameter £, the 
field strength through /3, and the extent of the dead zone 
through 9 W i n d. In addition, the mass flux f ma ss is used in 
the boundary condition of the poloidal momentum com- 
ponents. Boundary conditions at the stellar surface are 
identical as above, but now the dead zone has a zero mass 
flux, so that fmass{9). Note that in a completely analogu- 
ous way, we could allow for a latitudinal dependence of the 
stellar rotation rate fl*(9), or the magnetic field strength 

0(0)- _ 

This t = guess for an axisymmetric MHD wind 
is then time-advanced to a stationary solution. Figure 
shows the final stationary state, for the parameter val- 
ues Vesc — 3.3015c s », 7 = 1.13, a constant rotation rate 
corresponding to £ — 0.0156, (3 — 3.69, 6 W md = 60° 
and the mass flux in the wind zone set to the constant 
fmass = 0.01377, while it is zero in the dead zone. These 
parameters are as in the WD solution and the Saku- 
rai wind presented earlier. The initial field geometry has 
evolved to one where the open fieldlines are draped around 
a distinct bipolar 'dead' zone of limited radial extent and 
the prescribed latitudinal range. The outflow nicely traces 
the field geometry outside this dead zone. As seen from 
the figure, we have calculated the full poloidal halfplane 
and imposed symmetry boundary conditions at north and 
south pole. We used a polar grid of resolution 300 x 40 of 
radial extent [1, 50]r» with a radial grid accumulation at 
the base. The north-south symmetry of the final solution 
is a firm check of the numerics. The critical surfaces are 
also indicated in Fig. || and they differ significantly from 
the monopolar field solution shown in Fig. |i[ Again, at 
the polar regions, the Alfven and fast critical surface co- 
incides. Now, the B v also vanishes at the equator where 
conditions are such that the slow and the Alfven critical 
surfaces coincide. The B v component changes sign when 
going from north to south, as the rigid rotation shears the 
initial, purely poloidal bipolar magnetic field. This is dif- 
ferent from the Sakurai wind presented above, where the 
boundary condition on B v was taken symmetric about the 
equator. Note how the equatorial acceleration to super- 
Alfvenic velocities occurs very close to the end of the 
dead zone. The critical surfaces are all displaced inwards 
as compared to the monopolar case. 

Figure |B| shows that poloidal streamlines and fieldlines 
are parallel. The E v is below 3%. In Fig. ^, we show the 
latitudinal variation of the (scaled) density and the veloc- 
ity at two fixed radial distances in a polar plot. The space- 
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craft Ulysses and the on-board SWOOPS experiment pro- 
vided the solar community with detailed measurements of 
these quantities for the solar wind (McComas et al. 1998Q . 
Qualitatively, the measured poloidal density and veloc- 
ity variation resembles the curves from Fig. |^: the den- 
sity is higher about the ecliptic and there is a decrease 
in wind speed associated with the equatorial 'dead' zone. 
However, our computational domain extended to 50 stellar 
radii, while Ulysses measurements apply to larger radial 
distances. Note that we could use observed solar differen- 
tial rotation profiles, as well as mass fluxes and magnetic 
field strengths, to obtain a better MHD proxy of solar 
wind conditions. The extent of the solar coronal active 
region belt suggests the use of a 'dead' zone larger than 
modeled in Fig. |[ 

7. Conclusions and outlook 

We obtained polytropic stellar winds as steady-state tran- 
sonic outflows calculated with the Versatile Advection 
Code. We could relax an isothermal, spherically symmet- 
ric Parker wind, to a polytropic wind model. Subsequently, 
we included stellar rotation and a magnetic field, to arrive 
at the well-known Weber-Davis solution. We used fully 
implicit time stepping to converge to the steady-state so- 
lutions. The correctness of these ID wind solutions can be 
checked precisely. 

We generalized to 2D axisymmetric, unmagnetized and 
magnetized winds. Noteworthy is our prescription of the 
stellar boundary conditions in terms of the prescribed 
mass flux / mass and the way in which the parallelism of the 
flow and th e ficld lines in the poloidal plane is achieved. In 
Bogovalov ( 1996 ), the stellar boundary specified the nor- 
mal magnetic field component and the density at the sur- 
face, while keeping the velocity of the plasma on the stellar 
surface in the rotating frame constant. Our approach dif- 
fers markedly, since we impose the mass flux and ensure 
the correct rotational coupling of velocity and magnetic 
field. We refrain from fixing the density, as the analyti- 
cal treatment shows that the algebraic Bernoulli equation 
together with the cross-field momentum balance really de- 
termines the density profile and the magnetic flux func- 
tion concurrently, and should not be specified a priori. In 
fact, we let the density and all magnetic field components 
adjust freely at the base. This allows for the simultane- 
ous and self-consistent modeling of both open and closed 
fieldline regions, which is not i mmed iately possible when 
using the method of Sakurai ( 1985 ). By an appropriate 
initialization of the time-marching procedure used to get 
the steady-state solutions, we can find magnetized winds 
containing both a 'wind' and a 'dead' zone. 

The method lends itself to investigate thermally 
and/or magneto-centrifugally driven polytropic wind so- 
lutions. One could derive angular momentum loss rates 
use d in s tudies of s tellar rotational evolution (Keppens et 
al. 1995 , Keppens 1997 ). However, our immediate inter- 



N 




Fig. 5. Axisymmetric magnetized wind containing a 
'wind' and a 'dead' zone. Shown are the poloidal magnetic 
ficldlines and the poloidal flow field as vectors. Indicated 
are the three critical surfaces where Mf = 1 (dotted), 
M\ = 1 (solid line), and Mf = 1 (dashed). 



est is in the relaxation of the assumptions inherent in our 
approach. 

In this paper, we assume a polytropic equation of state 
throughout. All solutions are smooth and demonstrate a 
continuous acceleration from subslow outflow at the stel- 
lar surface to superfast outflow at large radial distances. 
Our polytropic assumption has to be relaxed to investi- 
gate the combined coronal heating/solar wind problem 
within an MHD context. This involves adding the energy 
equation. We plan to study possible discontinuous tran- 
sonic solutions containing shocks. We can then address 
the puzzling paradox recently raised by analytic investi- 
gations of translational symmetric and a xisym metric tran- 
sonic MHD flows (Goedbloed & Lifschitz |1997| , Lifschitz & 



10 



R. Keppens & J. P. Goedbloed: Numerical simulations of stellar winds 




al. (1996) are called for. Combined analytic and numerical 
studies of such axisymmctric steady-state flows have been 



South Pole 



0.5 



1.0 



Fig. 6. Polar plots of the scaled density (dashed) and ve- 
locity (solid) for two fixed radial distances: 10 and 20 
(thick lines) stellar radii. The 'dead' zone has a clear in- 
fluence on the latitudinal variation. 



Goedbloed |1997| , Goedbloed et al. |l998|) . The generalized 



Grad-Shafranov equation describing the cross-fieldlinc 
force balance has to be solved concurrently with the al- 
gebraic condition expressed by the Bernoulli equation. 
Rigorous analysis of the generalized mixed-type Grad- 
Shafranov partial differential equation, in combination 
with the algebraic Bernoulli equation, shows that only 
shocked solutions can be realized whenever a limiting line 
appears within the domain of hyperbolicity. Moreover, 



in Goedbloed & Lifschitz (1997) and Lifschitz & Goed- 
bloed (1997), it was pointed out that there are forbidden 
flow regimes for certain translationally symmetric, self- 
similar solutions of the MHD equations. The Alfven criti- 
cal point is in those solutions situated within a forbidden 
flow regime, which can only be crossed by shocks, ft is 
of vital importance to understand what ramifications this 
has on analytic and numerical studies of stellar winds, or 
on accretion-type flows where shocked solutions are rule 
rather than exception. Since the schemes used in VAC are 
shock-capturing, we have all ingredients needed to clarify 
this paradox. Numerical studies of self-s imilar solutions as 
those discussed in Trussoni et al. ( 1996 ) and Tsinganos et 



initiated in Goedbloed et al. (1998) and in Ustyugova et 
al. (^99|J). 

After those paradoxes are resolved, we will be in a 
position to relax the conditions of axisymmetry and sta- 
tionarity. While several authors have already initiated this 
daunting task (Gibson & Low|l998[ Guo & Wu [1998| Wu 



& Dryer 1997, Usmanov & Dryer 1995), we believe that an 
in-depth study of the subtleties involved with the various 
restrictions mentioned is still warranted. 
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